Astronomy & Astrophysics manuscript no. theory 


© ESO 2009 


July 14, 2009 





A theoretical framework for the detection of 
point-sources in Cosmic Microwave Background maps 

G\ ■ 
O 

! R. Vio 1 and P. Andreani 2 

<N . 

r—{ • Chip Computers Consulting s.r.L, Viale Don L. Sturzo 82, S.Liberale di Marcon, 30020 Venice, Italy 
K - 5 ' e-mail: robertovio@tin.it, 

i 2 ESO, Karl Schwarzschild strasse 2, 85748 Garching, Germany 
■ INAF-Osservatorio Astronomico di Trieste, via Tiepolo 11, 34143 Trieste, Italy 
e-mail: pandrean@eso . org 



^ : 

Or 



Received ; accepted 



O ' ABSTRACT 

Aims. The detection of point-sources in experimental microwave maps is a critical step in the analysis of the Cosmic Microwave 
i i Background (CMB) data. If not properly removed, these sources can have adverse effects on the estimation of the power- 
, spectrum and/or the test of Gaussianity of the CMB component. In the literature, various techniques have been presented to 
£N) ■ extract point sources from an observed image but no general consensus about their real performance and properties has been 
i> reached. Their characteristics have been studied essentially through numerical simulations based on semi-empirical models of 
V) , the CMB and the Galactic foreground. Such models often have different levels of sophistication and/or are based on different 

■ physical assumptions (e.g. the number of Galactic components and level of the noise). Moreover, the application of a given 
technique to a set of data (either simulated or experimental) requires the tuning of one or more parameters that unavoidably 
is a subjective operation. Hence, a reliable comparison is difficult. What is missing is a statistical analysis of the properties of 

' the proposed methodologies. This is the aim of the present paper. 

. Methods. The statistical properties of the detection techniques in the context of two different criteria, i.e. the Neyman-Pearson 

■ criterion and the maximization of the signal-to-noise ratio (SNR), are analyzed through an analytical approach. One-dimensional 
' as well two-dimensional signals are considered. The case of multiple observing frequencies is also addressed. 

Results. The conditions are fixed under which the techniques can work satisfactorily. Their limits are also illustrated and 
V~J ■ implementation details provided. We show that, exploiting some a priori information, it is possible to develop simple algorithms 
K% ' with performances not too far from those of more sophisticated but complex techniques. In this respect, a detection algorithm, 
tailored for high Galactic latitudes, is presented that could be useful in future ground-based experiments as, for example, the 
Atacama Large Millimeter/submillimeter Array (ALMA). 



Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 



1. Introduction 

The detection of point-sources embedded in a noise back- 
ground is a critical issue in the analysis of the experimen- 
tal Cosmic Microwave Background (CMB) maps. The es- 
timation of the power-spectrum of the CMB component 
and the test of its possible nonGaussian nature need the 
a priori detection and removal of these sources. In par- 
ticular the former operation is rather delicate. Given its 
importance, this subject has been extensiv ely considered 
in literature (see Herranz and Sanz~ll2008al and references 
therein). However, no widespread accepted conclusion has 
been reached. The reason is that most of the detection 
techniques presented in literature lack a sufficiently rigor- 
ous theoretical background. The statistical characteristics 



are derived from numerical experiments only. Especially 
in the experiments where it is necessary to simulate maps 
at different observing frequencies, this approach is haz- 
ardous for two reasons: 1) the models used by the various 
authors in the simulations are not the same; 2) when a 
detection technique has to be applied to a set of data (ei- 
ther synthetic or real) the tuning of some parameters is 
unavoidably subjective. Because of this contrasting results 
appear in literature. A safer procedure consists in study- 
ing the techniques in a well defined theoretical framework. 
Although some of the conditions assumed to fix the the- 
oretical context could be not realistic, anyhow it is pos- 
sible to obtain a more objective comparison as well as 
some useful indications that help to understand what it 
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is possible to expect when they are relaxed. The aim of 
the present paper is to provide a theoretical characteriza- 
tion of some of the detection techniques suited for both 
single and multiple-frequency CMB experiments. For ease 
of notation, initially the arguments will be developed for 
one-dimensional signals x = (x[Q],x[l], . . . ,x[N — 1]) T . 
Later they will extended to the two-dimensional situa- 
tion. The paper is organized as follows: in Sec. [2] the 
problem of point-source detection is described in detail 
for the case of a single spatial dimension. Here some 
standard material is presented with the aim to fix nota- 
tion and formalism. In Sec. [3] the situation is considered 
when more signals are available. In Sec.|4]a technique tai- 
lored for observation at high Galactic latitude is presented 
that will be useful in some experiments planned in the 
near future with innovative instruments as the Atacama 
Large Millimeter/submillimeter Array (ALMA). Finally, 
the conclusions are given in Sec. [5] In appendix [A] two 
procedures are presented that efficiently implement the 
techniques described in the text. These are extended to 
the two-dimensional case in appendix iBl 



2. Formalization of the problem 

The first step in the development of a detection technique 
is to fix the conditions that are pertinent to the problem 
of interest. In the case of CMB observations, the following 
conditions are commonly assumed: 

1. The point-sources have a known spatial profile s = ag. 
The amplitude "a" is a scalar quantity different from 
source to source, whereas g is a function, due to 
the instrument beam, that is identical for all of 
them. Function g is normalized in such a way that 
max{.g[0],.g[l],..., 5 [iV-l]} = l; 

2. The point-sources are embedded in a noise- 
background, i.e. the observed signal x is given 
by x = s + n. In other words, noise is additive. 
Usually, this is a reasonable assumption; 

3. Noise n is the realization of a stationary stochastic 
process with known covariance matrix 



equivalent to a decision problem where two hypotheses 
hold: 

Ho- x = n; 

Til x = n + s. ^ ' 

Under Hq the probability density function of x is given 
by p(x\Tto) whereas under Hi by p(x\TLi). At this point, 
it is necessary to fix the criterion to use for the detection. 
Clearly, one cannot hope to find all the sources present 
in a given signal. Hence, some choices are necessary. For 
example, one could decide that the non-detection or the 
misidentification of a bright source could be more impor- 
tant than that of a fainter one, or vice versa. A very com- 
mon and effective criterion is the Neyman-Pearson crite- 
rion that consists in the maximization of the probability of 
detection Pq under the constraint that the probability of 
false alarm Pfa (i-e., the probability of a false detection) 
does not exceed a fixed valu e a. The Neyman-Pearson the- 
orem (e.g., see iKav 1998 ) is a powerful tool that allows 
to design a decision process that pursues this aim: to max- 
imize Pd for a given Pfa = ol, decide Hi if the likelihood 
ratio (LR) 



p(x\Hi) 
p(x\H ) 

where the threshold 7 is found from 



Pfa 



{x:L(x)>~f} 



p(x\Ho)dx = a. 



(3) 



(4) 



C = E[nn T ]. 



(1) 



The test of the ratio ^ is called the likelihood ratio test 
(LRT). 

An important example of application of LRT is when 
noise n is Gaussian with correlation function C. Actually, 
in CMB experiments this condition is satisfied only for 
observations at high Galactic latitudes where the CMB 
emission and the instrumental noise are by far the domi- 
nant contributors. At lower latitudes, it is often assumed 
to hold locally. For example, the contribution to x of com- 
ponents that in small patch of sky presents linear spatial 
trends are often approximated with stationary Gaussian 
processes with a steep spectrum (e.g. 1/f noises). In any 
case, even if the Gaussianity condition was unrealistic, it is 
often made anyway since it allows an analytical treatment 
of the problem of interest and the results can be used as 
a benchmark in the analysis of more complex scenarios. 
With Gaussian n, it is 



Actually, because of the Galactic contribution, espe- 
cially at low Galactic latitudes, this hypothesis is not 
satisfied. However, it is assumed to hold locally. This 
allows the computation of statistics as the mean or 
the covariance matrix that, otherwise, should not be 
possible. Without loss of generality, it is assumed that 
E[n] = 0. 

Under these conditions, the detection problem consists in 
deciding whether x is a pure noise n (hypothesis Hq) or 
it contains also the contribution of a source s (hypothe- 
sis Hi). In other words, the source detection problem is 



with 



p(x\Ho) = Aexp 
p(x\Hi) = Aexp 

A = 



1 , 
— x 
2 



C~ L x 



— (x — s) 



^C-\x- 



(27r)Trdet5(C) 



The LRT is given by 

l(x) = ln[L(x)] = x T C^ 



s — 



"C -1 * > 7. 



(5) 
(6) 

(7) 

(8) 
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Hence, it results that Hi has to be chosen when for the 
statistic T(x) (called NP detector) it is 



T(x) = x 1 C~ l s > 7 




[a T C- l a\ 



= ot, 



i.e., 



Q- 1 (P^A.Ws T C- 1 8. 



(9) 



(10) 



(11) 



Here, Q(x) = 1 — $(a;) with Q(x) the standard Gaussian 
distribution function, Q _1 (.) is the corresponding inverse 
function. Equation (fTU|) is due to the fact that T(x) is 
a Gaussian random variable with variance s T C~ 1 s and 
expected values equal to zero under Tio and s T C~ 1 s under 
Til. For the same reason it is 



Pb = Q [Q- 1 (Pfa) - VsTC-'s) . (12) 
Equation (J9J can be written in the form 

T{x) = x T u > 7, (13) 

with 



u = C"V 



(14) 



From this equation appears that u can be thought as a 
linear filter of signal x, that is called matched filter (MF). 

2.1. Some comments on the use of the matched filter 
in practical applications 

There are some important points to stress with regard the 
MF when used in practical aplications. They are: 



T\{x) is a sufficient statistic ( Kav I Il998h . Loosely 
speaking, this means that T\(x) is able to summarize 
all the relevant information in the data concerning the 
decision ([2]). No other statistic can perform better. 
As a consequence, the claim that some filters (e.g. 
the Mexican hat wavelet, the scale-adaptive filters, 
the biparametric scale-adaptive filters...) are superior 
to MF according to the Neyman-Pearson criter ion 



(jBarreiro et al. I l2003t lLopez-Caniego et al. 1 120051 ) is 



not correct. It is the result of the use of im precise 
theoretical arguments (e.g. see I Vio et al. 1 120041 ): 



— It is worth noticing that if the amplitude "a" of the 
source is unknown, then Eq. ([9]) can be rewritten in 
the form 

T(x) = x T C 1 g > i, (15) 



with 7' = 7/a = Q^ 1 {Pfa)\J g T C~ x g . In other words, 
a statistic is obtained that is independent of "a", i.e. 
also in the case that the amplitude of the source is 
unknown, T(x) maximizes Pd for a fixed Pfa- The 
only consequence is that Pd cannot be evaluated in 
advance. In principle this quantity could be evaluated 



a posteriori by using the maximum likelihood estimate 
of the amplitude, a = x T C~ 1 g/g T C~ 1 g, but this is of 
little interest. More useful is that in real experiments 
one is typically interested in the detection of sources 
which have amplitudes characterized by a probability 
density function p(a). In this case, once Pfa is fixed 
to a value a and making to change "a" across the 
domain of p(a), the quantity 1 — Pd, with Pd as given 
by Eq. (fT2")l and s — ag, provides an estimate of the 
fraction of undetected sources as function of their 
amplitude; 

If s = Hs and x = Hx, then 

T{Hx) = x T C~ l s 

= x T H T H T C 1 H l Hs = T(x), (16) 

with H any invertible linear operator (matrix). A 
useful consequence of this property is that if signal 
x is convolved with a function (e.g., the beam of an 
instrument), this operation does not modify the opti- 
mality of MF. This fact could be useful in situations 
where more signals are available that are obtained 
with different point spread functions (see below); 

The arguments above are developed under the implicit 
assumption that sources do not overlap and that their 
position is known in advance. In practical applications 
the first condition is satisfied - strictly speaking - 
only in very high resolution maps and it is assumed to 
be always valid for those sources above the confusion 
noise. The second statement, i.e. that the position of 
the sources is known in advance, is not true and the 
standard procedure consists in filtering x by means 
of u and in computing T(x) for the peaks in the 
resulting signal. Although there is no guarantee that 
a peak in the filtered signal marks the true position 
of a source even in the case this is effectively present, 
theoretical arguments as well as years of application 
in real-life problems have proved that this procedure 
is rather robust and able to provide excellent results; 

If the Gaussianity of n is relaxed, then MF is no longer 
optimal in the Neyman-Pearson sense. However, it re- 
mains optimal with respect to the signal-to-noise ratio 
(SNR)li This means that, independently of the nature 
of the noise, MF provides the greatest amplification of 
the signal with respect to the noise. This can be easily 
verified through the minimization of the variance of the 
filtered noise u T n with the constraint that u T s = a 
(i.e. filter u docs not modifies the amplitude of the 
source), in formula d 

m snr = argmin[u, T Cit — \(s T u — a)] (17) 



1 Here, the quantity SNR is defined as the ratio between the 
squared amplitude of the filtered source with the variance of 
the filtered noise. 

2 We recall that the functions "arg min F(x)" and 
"argmax_F(a;)" provide the values of x of for which the 
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with A a Lagrange multiplier. It is not difficult to see 
that 

msnr = aC- x s/[s T C- x s], (18) 
= C 1 g/[g T C 1 g], (19) 

i.e. apart from a normalizing factor, msnr is given 
by Eq. (|14p . Since msnr is optimal as concerns the 
maximization of the SNR, no other methods can out- 
perform it in this respect. For this reason, expedients 
as the introduction of additional constraints and/or 



of free parameters i n u (e.g. see ISanz et al. 1 12001 



Herranz et al. I 120021 ) has the only effect to reduc e 



the detection performances (e.g. see lVio et al. 2002 ). 
One of the benefit in using «snr is that the value of 
u|nr s provides directly an unbiased estimate of the 
amplitude "a" of the source. However, it is necessary 
to stress that in practical applications, where the true 
position of the source is not known and it is necessary 
to apply the procedure described above, this is no 
longer true; 

When s is a long signal, some computational prob- 
lems come out. In fact, if the size of s is N, then C 
is a N x N matrix. Hence, the computation of the 
quantity x T C~ 1 s can become quite expensive. The 
computational burden can be alleviated if C, that is a 
Toeplitz matrix, is approximated with a circulant ma- 
trix C. This is because C can be diagonalized through 



(N)CFfN) 



(20) 



Here, F(n) is the N x N one-dimensional Fourier ma- 
trix that is a complex, unitary, and symmetric matrix 
whose elements are given by 



( F (N))k 



kl 



-27Tt(fc-l)(i-l)/JV 



N 



(21) 



denotes the one-dimensional Fourier 



= \J — 1, F N is the complex conjugate 



Symbol "~" 
transform, i 

transpose of Ftf, and S is a diagonal matrix contain- 
ing the eigenvalues of C (i.e., the power-spectrum of 
n). After that, since Fn is a unitary matrix, one ob- 
tains that 



x T C- l s « (x' F? N) )(F (N) CF? N) y\F iN)S ) = 



(22) 



= x H ~E s = x"{BlAG[E "]©5). (23) 

Symbol "0" denotes the element-wise multiplication, 
and DIAG[Z] a column vector containing the diagonal 
elements of the square matrix Z. In obtaining this 
result, the fact that F^C^Ff^ = (F^CFf^)- 1 



-H 



has been used. In principle, all the terms in Eq. (|23l) 
could be computed quite efficiently by means of 
fast Fourier transform (FFT). In particular, array 



function F(x) has the smallest and greatest value, respec- 
tively. 



DIAG[S ] can be obtained by means of the re- 
ciprocal of the FFT of the autocorrelation function 
c(t) = E{n[fc + r]n[fc]}. Actually, in using C, there is 
the problem that both s and x are implicitly assumed 
to be periodic signals with period N. In general, this is 
not true. Hence, boundary effects are to the expected 
in the computation of the FFT. However, since s 
typically has finite spatial support, these effects can 
be easily avoided by padding it with a sequence of 
leading zeros longer than the correlation length of the 
noise. This is visible in the last three panels in Fig. [T] 
where the array C _1 s, that approximates the MF 
u = C~ 1 s in Eq. (fT4| . is shown when C is constructed 
with the correlation function shown in the first panel 
and s is a rectangular function (shown in the same 
panel) that is padded with an increasing number of 
leading zeros. When this number is sufficiently large, 
it is evident that C~ 1 s « C~ l s. In principle, the same 
method could be applied to x. However, because of the 
noise, usually this signal has no finite spatial support. 
Hence, in order to evaluate x T C~ 1 s, often one wishes 

to compute the inverse FFT of DIAG[£ ] s, 
to remove a number of leading elements from the 
resulting array equal to that of the zeros used in the 
padding operation, and then to calculate the scalar 
product with x. 

In the case n is a colored noise (i.e. C is not a diagonal 
matrix), then the length N oi s and x should be longer 
than the correlation length of n. This is clearly visible 
in Fig. [3 where the different performances of MF are 
compared when n is the realization of a Gaussian pro- 
cess with a correlation length of about 100 pixels and s 
is a Gaussian with a = 1 , dispersion set to three pixels 
that is computed, respectively, on 13, 101 and 301 pix- 
els. It is evident that when N > 100 the performance 
of MF becomes independent of this parameter. The 
comparison is based on the so called receiver operating 
characteristics (ROC) that is a plot of Pd versus Pfa- 
This kind of plot is very useful since it allows a direct 
visualization of the detection performances of a given 
technique. More specifically, the ROC should always 
be well above a 45° straight lineline since this corre- 
sponds to a detection performance identical to that of 
flipping a coin, ignoring all the data. 



3. Extension to the mutiple-frequency case 

In the context of CMB observations, there is a further 
complication in that there are M signals x^ = + n^, 
&k = a kQki k = 1,2, ...,M, coming from the same sky 
area that are taken at different observing frequencies. 
Here, ah is the amplitude of the source at the fcth ob- 
serving frequency, whereas g k is the corresponding spatial 
profile. For ease of notation, all the signals are assumed to 
have the same length N. In general, the amplitudes {a^} 
as well as the profiles {g k } are different for different k. 
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X = 




T 

x 2 i ' • 


! x mJ 


s = 


1*1, 


s 2 ! • • • 


,T iT 

> *mJ i 


n = 




T 

n 2 , . . 


T iT 
■ I n M\ 



However, if one sets 

(24) 
(25) 
(26) 

it is possible to obtain a problem that is formally identical 
to that treated in the previous section. Hence, the MF is 
still given by Eqs. (|13 p - (|14p and is named multiple matched 
filter (MMF). The only difference with the classic MF is 
that now C is a (NM) x (NM) block matrix with Toeplitz 
blocks (BTB): 



The solution is 

«snr = C-'G^C-'Gy 



(32) 



It is evident that, contrary to MMF, with msnr it is pos 
sible to obtain a statistic Tsnr(:e), 

Tsnr(x) = s t m S nr, 



(33) 



that is independent of the unknown source amplitude a. 
The price to pay is a reduced detection capability. In fact, 
better results should be obtainable if the SNR of each 
signal was maximized. However, this is an operation that 
requires the knowledge of a. For msnr, it is 



Ci 



M 



(27) 



Pfa = Q 



7 



[1 T (G T C~ 1 G)- 1 1] 1 / 2 



a, 



(34) 



C 



Ml 



• MM 



i.e. each of the Cy blocks is constituted by a N x N 
Toeplitz matrix. In particular, Cu — E[n^n^] provides , 
the autocovariance matrix of the ith noise, whereas = Pd = Q [ Q (Pfa) 
E[n»nJ], i ^ j, the cross-covariance matrix between the 
ith and the jth ones. 

In spite of these similarities, when M > 1 some difficul- 
ties arise. In particular, T(x) cannot be written in a form 
equivalent to Eq. (fT5|) . This has important consequences 
in the fact that if the amplitudes {a^} are unknown, then 
T(x) cannot be computed. In other words, if the spectral 
characteristics of the radiation emitted by a source are 
not fixed, then the MMF is not applicable. This forces to 
resort to an approach based on the maximization of the 
total SNR of the filtered signals. Following this approach, 
model (fl~7|) has to be modified in the form 



that again is a quantity independent of the source ampli- 
tude, whereas as expected the same is not true for Pd 



i t 1 



[l T (G T C _1 G)- 1 l] 1 /2 



(35) 



In two recent work s iHerranz and Sanz I (|2008af ) and 
Herranz et al. I (|2008bl) have proposed a "new" class of 



filters, the so called matrix filters. Their idea is to filter 
separately each of the signals xy. in such a way to ob- 
tain unbiased estimates of {a^} and at the same time to 
simultaneously minimize the total variance of the filtered 
signals. In the spatial domain, this problem can be written 
in the form 

C/snr = argmin{Tr[C7 T GJ7 - A(G T U - I)}}, (36) 
u 



msnr = argmin[it T G« — \ T (S T u — a)], (28) 



where u = [u± , u 2 , 
A = [Ai, A2, 
(NM) x M matrix 



Am] T , a 



S = 



uJ f ] T is an (NM) x 1 array, 
= [ai, d2, ■ ■ ■ , aj\,f] T and S is a 



(29) 



with "Tr" the Trace operator, 

( Mil 
"21 



u 



U12 
U22 



UiM \ 

U 2 M 



(37) 



( S1 


. 


. 


\ 





«2 


. 








. 




V 


. 


■ s M 


) 



\ UM1 U M 2 ■ ■ ■ U M M ) 

a (NM) x M matrix of filters and 



A = 



/ An A12 
A21 A22 



Aim \ 

A2M 



(38) 



with = [0,0,..., 0] T a N x 1 array. Now, since a = 
diag[a]l and S T = diag[a]G T with 1 = [1, 1, ... , 1] T and 
diag[a] a diagonal matrix whose diagonal contains a, it is 
trivial to show that Eq. (|28|) is equivalent to 



msnr — argmin[M T Gtt — A, (G u—1)], (30) 

u 

where A* = (diag[a]) _1 A and 
(9i 

G= " U1 ' - ° . (31) 



g 2 
\ 



. 

• 9m J 



\ Ami Am2 • ■ ■ Amm / 

a M x M matrix of Lagrangian multipliers. The solution 
is 

t/sNR = C- l G[G T C- l G]-\ (39) 
Through Csnr it is possible to define a set of M statistics 

T Ma t F (a;) =a; T [/ SN R (40) 

that can be considered individually. However, the fact that 

msnr = U snrI (41) 

indicates that msnr and U snr essentially represent the 
same filter. The only difference is that, after filtering sig- 
nals {xk}, the latter does not compose them together 
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as the former does. We call «snr the modified multiple 
matched filter (MMMF). 

In appendix [X] two procedures are presented that work 
in the Fourier domain and that allow a very efficient com- 
putation of both msnr and U snr- In appendix [Bl the 
methods are extend to the two-dimensional case. 

3.1. Alternative techniques 

A strategy to deal with multiple-frequency images that is 
alternative to the approaches presented above consists in 
composing signals {xk} together in a single array y. After 
that, the classic MF could be applied. Without a priori 
information, the most obvious method is 



M 



y 



J2 H 



(42) 



where Hk is an operator (matrix) such as HkQk = 9 inde- 
pendent of "fc" . In this way, the effects due to the different 
instrumental beams can be avoided. As seen above, the use 
of Hk does not represent a problem since it does not mod- 
ify the statistic T(x). This approach, that we indicate as 
summed-image matched filter (SMF), can be expected to 
provide results close to those of the MMF when the am- 
plitudes {afe} as well as the level of the noises are similar. 
However, if this condition is not satisfied, its performances 
rapidly worsens especially if rik and n;, k ^ I, have some 
degree of correlation. This last condition is typical of CMB 
signals where rik = n c + ek with n c a component inde- 
pendent of k (i.e. the CMB contribution), and e& a white- 
noise process due to the electronic of the instrument. In 
this case, it could be preferable a weighted sum as 



M 



y = y^WkH k Xk, 



fe=l 



where 



w T l = 0, 
w T w = 1. 



(43) 



(44) 
(45) 



with w = [w\, W2, ■ ■ ■ , Wk] T ■ The first constraint implies 
that the contribution of n c in x is completely removed, 
whereas the second one provides a normalizing factor. We 
indicate this method as weighted matched filter (WMF) 
and the particular case where w = [p,p,..., —(M — l)p] 
with p = 1/ y(l — M) + (1 — M) 2 , as uniformly weighted 
matched filter (UWMF). This last corresponds to a situa- 
tion where only one signal is used to eliminate the shared 
noise component n c , whereas the others are given an iden- 
tical weight. The UWMF is not very effective, but it can 
be useful in absence of a priori information on the spectral 
properties of the sources. 

The performance of these methods depends critically 
by many factors as, for example, the relative value of the 
amplitudes {a^}, the correlation lengths of the rel- 

ative importance of the noises {e^} as well as the degree 



of correlation between the different observing frequencies. 
This fact is evident in Fig. [3] that shows the Pfj vs. the 
amplitude a 2 when M = 2 and Pfa = 0.01 and 0.1, re- 
spectively. Here a% = 1, g is a Gaussian with dispersion set 
to three pixels, rife = n c + with n c a zero-mean, unit- 
variance Gaussian process whose autocorrelation function 
has Gaussian profile and dispersion set to ten pixels and 
finally ek, k = 1,2, two independent Gaussian white-noise 
processes. Two different cases are considered for the noises 
efc. In the first the noises e± and e 2 have the same disper- 
sion, i.e. <7i = (72 = 1, whereas in the second one o\ = 1 
and (72 = 0.5. When M — 2, the only possible weights 
for the WMF are e ither u = [l/y^; -l/VW . This is the 
case considered by IChen and Wright ] (|2008h . 

One indication that comes out from these examples is 
that, as expected, MMF outperforms all the other meth- 
ods. Moreover, unless the level of noises {ek} are similar, 
MMMF outperforms SMF. Heuristically, these results can 
be explained by the fact that through MMF a sum of 
signals is computed that is weighted by means of both 
the intensities {a^} and the noise levels {o~k}, whereas 
with MMMF the weighting is based on the noise levels 
only. In SMF there is no weighting at all. A different sit- 
uation is for WMF. Here the two signals are subtracted. 
This operation provides a signal for which i) the correlated 
part in the noises is zeroed; ii) the total amplitude is 
a = a\ — a2', hi) the variance a 2 of the instrumental noise 
is a 2 = of + of. Since SNR = a 2 /a 2 , it is not difficult to 
realize that a benefit in the detection capability happens 
when a-i <C a\. In the case case M > 2, similar arguments 
indicate that good results can be expected when a channel 
"Z" is available with a/ <C a^i ■ This is the case expected 
in CMB applications (see below). 



4. An approach for CMB experiments at high 
Galactic latidude 

In the near future some innovative ground-based ex- 
periments are planned for very high spatial resolution 
observations as, for example, with the Atacama Large 
Millimeter/submillimcter Array (ALMA). One important 
advantage of these experiments is that, contrary to the 
satellite observations, they will allow a certain control of 
the experimental conditions. Hence, a full exploitation of 
the capabilities of this facility (and other instruments) re- 
quires a careful planning of the observations. The record- 
ing of good quality data will allow an effective application 
of the chosen techniques for the detection of point-sources. 
For instance, with instruments as ALMA it is possible to 
plan observations at high latitude fields to map sources 
at very high spatial resolution in sky regions dedicated to 
CMB observations. In this case a detection technique is 
necessary that takes into account such a specific experi- 
mental situation. 

In the context of point-source detection, data can be 
thought as two-dimensional discrete maps each 
of them containing N p pixels, corresponding to M different 
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observing frequencies (channels), with the form 

X i =S i +M i . (46) 

Here, Si correspond to the contribution of the point- 
sources at the ith frequency, whereas Mi denotes the cor- 
responding noise component. At high Galactic latitudes, 
the CMB component is expected to be the dominant one. 
Hence, Mi may be modeled with 



Mi = B + £i 



(47) 



where £i is the experimental noise corresponding to the 
ith channel and B the contribution of the CMB compo- 
nent that is frequency-independent. The contribution of 
the point-sources is assumed in the form 



a,iG, 



(48) 



with ai the amplitude of the source to the ith channel. 
According to Eq. (j48| . and without loss of generality, all 
the sources are assumed to have the same profile Q inde- 
pendently of the observing frequency. In fact, although in 
general this will not be true, as written in Sect. 12. li it is 
possible to meet this condition by convolving the images 
with an appropriate kernel with no consequences. In the 
following it is assumed that the components {Si} are the 
realization of stationary, zero-mean, stochastic processes. 

The main feature of model (gS])-(|17]) is that the CMB 
contribution does not change with frequency. Hence, fol- 
lowing the suggestion in Sec. 13. 1\ the WMF can be used. 
Of course, in order this approach be effective, it is neces- 
sary that the amplitude of a given point-source is not the 
same in all the observing channels. 

In the case M = 2, (i.e. only two maps are avail- 
able), the only possible solution is w T = [l/\/2; —1/^/2] or 
[-1/V2; 1/V2]. However, for M > 2 more degrees of free- 
dom are available. This allows the selection of the weights 
in such a way that specific conditions are satisfied. In par- 
ticular, one could wish that, after the linear composition 
of the maps, the quantity 



R(w\a) 



(wTDw) 1 / 2 ' 



is maximized, i.e. 



w = argmaxi?,(«;|a). 



(49) 



(50) 



Here, D is the Mx M cross-covariance matrix of the noise 
processes whose (i,j)th entry (-D)y is given by (D)ij = 
E[VEC T [£ l ]VEC[£ J ]]/N p , with VEC[5] the operator that 
transforms a matrix £ into a column array by stacking its 
columns one underneath the other. The rationale behind 
this choice is the the quantity R(w\a) is a measure of the 
amplitude of the point-source in the weighted map with re- 
spect to the standard deviation of the measurement noise. 
The larger R(w\a) the more prominent is the point-source 
with respect to the noise background; an attractive situ- 
ation in problems of source detection. The maximization, 
via the Lagrange multipliers method, of R(w\a) with the 



constraint (|44|) provides the following system of non-linear 
equations 

{I ^ll T )(aw T Dw Dwa T w) = 0. (51) 

It can be solved through an iterative algorithm based on 
the Newton method 



w k+ i 



w k + Aw k 



where 
and 

with 



Ji(w) 



\w k + Aw k \ 7 
Aw k = -J~ 1 (w k )R(w k ) 

4 

j{w k ) = jj(w k ), 

2=1 

2aw T D; 

—a T wD — Dwa T ; 

2 



M 



ll 1 aw 1 D\ 



Ji(w) = 



^-(a T wll T D + 11 T Dwa T ). 



(52) 
(53) 
(54) 

(55) 
(56) 

(57) 
(58) 



Here, three points are of concern: a) the iteration can be 
initialized with a starting guess Wq that contains random 
entries satisfying the constrains (j44"|) - ([4"5|) ; b) The normal- 
ization term in Eq. (|52[) implements the constraint (|45|) . 
This cannot be done via a Lagrange multiplier since the 
array w that maximizes the quantity R(w\a) can be 
determined unless a multiplicative constant. Hence, the 
Lagrange multiplier corresponding to the constraint (I45| 
should be equal to zero; c) For the same reason, matrix 
J(w k ) is rank deficient, hence J~ 1 (w k ) has to be under- 
stood as Moore-Penrose pseudo-inverse. 

Once that the map y has been produced, the MF can 
be used with no necessity to take into account the char- 
acteristics of the CMB. This is particularly useful in sit- 
uations where only small patches of sky are available and 
hence the sizes of X are much shorter than the correlation 
length of M . As indicated earlier, we call this method the 
weighted matched filter (WMF). 

As for MMF, the procedure described above needs that 
the array a be specified. In fact, if the weights w are com- 
puted for a source with a given a, they will be optimal only 
for any other source with amplitude oc a. However, it is the 
general trend that does matter. For example, satisfactory 
results can be expected for the sources that present steep 
spectra with similar behaviors (see below). This means 
that the optimization of w can be carried out for subsets 
of sources for which the above condition approximately 
applies. 

4.1. Numerical experiments 

In this section we present some numerical experiments to 
test the performances of WMF. In particular, we consider 
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a scenario where three different observing frequencies are 
available. We make the simplifying assumption that all 
the channels have the same point-spread function (PSF) 
which is a two-dimensional circular symmetric Gaussian 
normalized to have a peak value equal to one and with a 
dispersion set to three pixels. In correspondence to the ith 
channel the contribution Si due to a point-source is Si = 
ciiG, whereas the terms C„ in the covariance matrices C 
are given by 

C u = C B + <y\l. (59) 

Here of is the variance of the instrumental noise, assumed 
of Gaussian white-noise type, and C b is the covariance 
matrix of the CMB component sampled with a step of 
3.52' on a regular two-dimensional grid. For i ^ j, it is 
Cij = C b- This scenarios mimics that expected for the 
" Low- Fr equency Instrum ent" mounted on the PLANCK 
satellite (jVio et al. 1120031 ) . The available data are assumed 
in form of square maps containing (101 x 101) pixels each 
0. This size is large enough, with respect to the correla- 
tion length of C, to make results independent of it. For a 
power law spectrum as S=av a the amplitude of a point- 
source at an observed frequency 1/2, given its amplitude at 
an observed frequency V\, expressed in Thermodynamic 
temperature can be written as: 



where 



(3 = 



(el 3 - l) 2 

hv 
~KT 1 



(60) 



(61) 



(62) 



with h the Planck constant, K the Boltzmann constant, 
T the CMB temperature, and a is the spectral index. 
When expressed in antenna temperature this latter takes 
the value 1.6 for the infrared sources and —3 for the ra- 
dio ones (dominated by synchrotron emission). Note that 
this is an approximated expression, usually used in CMB 
experiment to make a direct comparison between the sky 
temperature brightness and the sources brightness. 
Fig. g] compares the ROC for MMF and WMF for, respec- 
tively, the radio and the infrared sources case, in a situa- 
tion of relatively low SNR and with the level of noise that 
is the same for all the channels. The amplitudes <2j, com- 
puted through Eq. (|60jl. correspond to the 30, 44, 70 GHz 
observing frequencies. We allow a 10% deviation from the 
amplitudes given in (f6T)|) . Hence, in the same figures the 
90% confidence envelopes are shown for both MMF and 
WMF that have been obtained by computing the ROC for 
a set of one hundred arrays a, = a+Aa iy i = 1,2, . . . , 100, 



3 In practical CMB applications, the fact to working with 
square patches of sky is not a limit since, independently of the 
shape of the available maps, point-source detection is typically 
carried out on small spatial windows sliding across the sky. 
This is for computational reasons as well as because the noise 
contaminating the maps has no uniform spatial characteristics. 



Adi a Gaussian random array with mean zero and co- 
variance matrix 0.1a T J. For reference, the result obtain- 
able with the UWMF is also plotted. The improved per- 
formances of WMF with respect to UWMF is evident. 
Moreover, although as expected the MMF is always supe- 
rior to WMF, their performances are rather similar. The 
reason is that, in the case of sources with steep spectra 
as given by Eq. (|60|) , there are at least two observing fre- 
quencies, say "Z" and "fc", for which a; <fli. As seen at 
the end of Sec. [3l this is the condition for WMF to work 
well. 

5. Conclusions 

In this work, the problem of point-source detection in noise 
background has been addressed from a theoretical per- 
spective. This allowed an objective comparison of the ex- 
pected performance of the various techniques. From this 
comparison it is evident that "in se" no method is su- 
perior to the other ones. The differences are due to the 
amount of a priori information that they exploit. In other 
words, it appears that the effectiveness of a technique is 
not linked to its sophistication rather to the ability in us- 
ing the available a priori information. In particular, the 
methods based on the Neyman-Pearson criterion can be 
expected to provide better performances than those based 
on the maximization of the signal-to-noise ratio because 
they make use of the amplitude of the sources. For the 
same reason, an improvement in the detection capability 
can be expected when more signals are available that cor- 
respond to the same sky area taking at different observing 
frequencies. On the other hand, the a priori information 
that at high Galactic latitudes the dominant components 
are the electronic noise and the CMB, with this last inde- 
pendent of the observing frequency, suggests that a subop- 
timal approach based on an opportune linear combination 
of the signals that eliminates the CMB contribution could 
have a performance close to MF. Hence, it is useful in 
practical applications as it avoids the estimation of the 
covariance matrix (or the power-spectrum) of the CMB. 

As last remark we would like to add the following. The 
superiority of a theoretical approach does not diminish 
the usefulness of the numerical experiments. It is, how- 
ever, a bad habit to fix the characteristics of a statistical 
methodology only by means of numerical simulations. The 
application of a detection technique to a set of (either real 
or synthetic) data requires the tuning of some parameters 
that is a subjective operation. As shown in Sec. 12.11 MF 
requires that the position of the candidate source is known 
in advance, while in practical application this piece of in- 
formation lacks. The common solution consists to filter x 
with u and then to compute the statistics T(x) for the 
peaks in the resulting signal, but, because of noise, there 
is no guarantee that a peak in the filtered signal marks 
the true position of a source even in the case this is ef- 
fectively present. In a numerical experiment this implies 
the definition of a window, around the true position of 
the source, where a peak is assumed to identify a source 
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candidate. The size of the window is a parameter that can 
have important consequences but there is not an objec- 
tive criterion to fix it. The same holds also for the other 
techniques. Moreover, as remarked again in Sec. 12.11 (see 
also Fig. [2]) , the performance of MF depends on the rela- 
tive size of the spatial region where the signal is sampled 
with respect to the correlation length of the noise. Hence 
the comparison of MF with other filters can give different 
results according to size of the patch of sky that are consid- 
ered. For these reasons it is risky the use of classes of filters 



as the mexican hat wavelet family (jLopez-Caniego et al. 



l2006a| bl) that lack any theoretical justification and whose 
effectiveness is supported only by numerical experiments. 
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Fig. 1. Experiment that shows how the matched filter u = C~ 1 s, Eq. (|14[) . with C the covariance matrix defined in 
Eq. {T]) (a Toeplitz matrix), can be well approximated using a circulant matrix C if signal s is padded with a number 
of leading zeros sufficiently large (see Sec. I2.ip . Top-left panel: signal s (a rectangular function) padded with (blue 
line), 20 (red line) and 100 (green line) leading zeros. For reference, the covariance function (cyan line) is shown that 
is used to form the Toeplitz matrix C and its circulant approximation C. The top-right, bottom- left and bottom-right 
panels compare array C~ 1 s with C~ 1 s for the three zero-padding situations mentioned above (respectively, 0, 20 and 
100 zeros). It is evident that when the number of zeros is larger than the correlation length of the noise n then C _1 s 
and C 1 s are almost indistinguishable. 
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Fig. 2. Receiver operating characteristics (ROC) of 
Matched Filter (MF) for a noise n with a correlation legth 
of about 100 pixels and three different lengths of signals s 
and x (see text end of Sec l2.1|) : blue, red and green lines 
correspond to 13, 101 and 301 pixels, respectively. The 
45° straight line represents a poor detection performance 
which is identical to that of flipping a coin, ignoring all 
the data, i.e. better performances correspond to lines well 
apart from this one. It is evident the bad performance of 
MF when signals are used that are shorter than the cor- 
relation length of noise. 
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a, = 1; o, = 1; o 2 = 1; P FA = 0.01 a, = 1; o,=1; o 2 = 1; P FA = 0.1 
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S a 2 

Fig. 3. Comparison of the performance of the multiple-frequency detection techniques described in Sec. [3] in the case 
of two observing frequencies: Probability of detection (i"b) vs. the amplitude 02 of the source signal at frequency v-i 
when the amplitude of the first one at frequency v\ is a\ = 1. For the right panels the probability of false alarm (Pfa), 
i.e. the probability of a false detection, is fixed to 0.01, whereas for the left ones is fixed to 0.1. For the top and the 
bottom panels, the standard deviation of the noise in the two signals is set to a% = a% = 1 and a% = 1 and &i = 0.5, 
respectively. Here, MMF = multiple matched filter, SMF = summed-image matched filter, WMF = weighted matched 
filter, MMMF = modified multiple matched filter. NB. In the top panels MMMF and SMF are perfectly overlapping. As 
expected, the superiority of MMF is unquestionable (it provides the best theoretical performance). The performance 
of the other filters depends on the relative importance of a-i with respect to a\. 
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Fig. 4. Results from the numerical experiment described in Sec. I4.ll relative to infrared sources (IRS) (left panel) 
and radio sources (RS) (righ panel). Comparison of the ROC of the multiple matched filter (MMF) with those of the 
weighted matched filter (WMF) and of the uniformly weighted matched filter (UWMF). MMF is used as benchmark 
since it provides the best theoretical detection performance. Here, ai and Gi (i = 1 — > 30 GHz, i = 2 — > 70 GHz, 
i = 3 — > 100 GHz) are respectively the amplitude of the source and the standard deviation of the instrumental 
Gaussian white-noise corresponding to the ith observing frequency in units of the standard deviation of the CMB 
signal. The 90% confidence envelopes are shown for both MMF (gray, dot-dashed lines) and WMF (green lines) that 
have been obtained by computing the ROC for a set of one hundred arrays a + Aa, with a = [1.00; 0.50; 0.11] for 
the IRS and a = [1.00; 0.33; 0.03] for the RS, Aa a Gaussian random array with mean zero and covariance matrix 
C = Q.la T I. For WMF the plotted ROC corresponds to the weights w w [0.74; -0.06; -0.68] T for the IRS and 
w = [0.78; -0.18; -0.60] T for the RS. It is evident that MMF and WMF provide similar performances. 
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Appendix A: Efficient numerical implementations 
of the MMMF: one dimensional case 

Filters usnrOe) an d Usnr(x) given, respectively, by 
Eqs. (|32| and Eq. (|39"|) can be efficiently computed in the 
Fourier domain. If in the BTB matrix C given by Eq. (|2T|) 
the Toeplitz blocks are approximated with circulant ma- 
trices Cij , then each of them can be diagonalized by means 
of the (NM) x (NM) matrix 



then 



5 = I( M ) <S> F {N) , 

with "(g)" the Kronecker product and I(m) 
identity matrix. In this way a matrix 



El 



£ = 



-'JVl 



A I 



J MM 



(A.l) 

the M x M 

(A.2) 



is obtained where are diagonal blocks containing the 

eigenvalues of Cy. The quantities G 1 C^ 1 and G T C~ 1 G, 

can be computed firstly by solving the system of linear 

~ ~h ~ 

equations SZ = G, and then through the product G Z 

(or through the procedure described in Sec. l2.ip . Since ma- 
trix £ is highly structured and sparse, the computational 
load is not excessive. This approach is suited to be imple- 
mented in high-level programming languages as MATLAB 
that allow a friendly handling of arrays and matrices. 

A more efficient algorithm, but suited for low-level pro- 
gramming languages as C or FORTRAN, is obtainable re- 
arranging the elements of a;, s, n and u according to the 
so called row rollout order, i.e., 

x=\x[0],x[l},...,x[N-l]} T , (A.3) 

with x[i] = [xi[i], X2 [*],••• , £Af [*]] T , and similarly for s, n 
and u. After that, models. ([30]) and ([36]) . as well as the 
corresponding solutions (|3"!?|) and (|3"9")l . still holds with C 
and G replaced, respectively, by 



/ Q[0] C[-l] ... C[-(N-1)]\ 
C[l] C[0] ...C[-(N-2)] 



(A.4) 



\ C[N - 1] C[N - 2] ... C[0] / 
where C[t] — E[n[i]n T [i + r]], and 

G=[DCS [£ 1 ],DCSi[g 2 ],...,DCS M _ 1 [gJ], (A.5) 
with 



£ fe = [5 fc [0],Of M _ 1) ,3 fc [l],Of M _ 1) ,...,g fe [7V-l],Of M _ 1) ] T . 

(A.6) 

Here, DCS; [.] denotes the down circulant shifting operator 
that circularly down shifts the elements of a co lumn array 
by I positions. Now, it can be shown again (see|K av II 19981 
page 504) that, if one sets 



(A.7) 



where S is a block diagonal matrix 

/So o 

s = 



(A. 



(A.9) 



with 



Vo 

/ Puifi) PM) 

\P M l{fi) PM2{fi) 



iJV-1 



PlM(fi) 



(A.10) 



MM 



(fi) 



fi = i/N, i = 0, 1, . . . , N— 1. Here, Pu(fi) represents the 
cross-power-spectrum at frequency fi between and n; . 
Similarly, for a given signal r , the array 5^11 provides the 
corresponding FFT, 

/ m \ 

£[/i] 

^r = r= . (A.ll) 

From these considerations, it results that 

^snr = ^QiG 1 ^ X G)-\ (A.12) 

with 

G = [DCSo^LDCSigy, . . . ,DCS M -i[9 M ]]. (A.13) 
Here, the advantage is represented by the fact that 

/r o \ 



V o 



(A.14) 



iN-l 



i.e., the inversion of S, that is a large matrix with size 
(NM) x (NM), can be obtained through the inversion 
of a number TV of much smaller M x M blocks. This 
fact, coupled with the structure of G, allows to com- 
pute efficiently IZsnr through block-matrix operations. 
Equation (IA.12|) provides t h e discr ete ve rsion of the re- 
sult by iHerranz and Sariz" ( 2008a ) and Herranz et al. 



(|2008bf ). Finally, as for the spatial domain, it is 

"SNR^IZsNR 1 - (A.15) 



Appendix B: Extension of MF and MMMF to the 
two-dimensional case 

B.l. Single-frequency observations 

The extension of MF to the two-dimensional signals X 
and <S is conceptually trivial. If one sets 



s = VEC[5]; 
x = VEC[Af]; 
n = VEC[Af], 



(B.l) 
(B.2) 
(B.3) 
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formally the same problem is obtained as that given by 
Eq. However, again some computational issues come 
out. The question is that, even for moderately sized sig- 
nals, matrix C becomes rapidly huge. In fact, if S con- 
tains N p pixels, then C is a N p x N p matrix. Similarly to 
the one-dimensional case, the computational burden can 
be alleviated by resorting to the Fourier domain. If <S is a 
N r x N c rectangular map, then C is a matrix of type block 
Toeplitz with Toeplitz blocks (BTTB), and it can approxi- 
mated with a matrix of type block circulant with circulant 
blocks (BCCB). In fact, a BCCB matrix C can be diago- 
nalized through 

(B.4) 



where T = F^ r ^F^ c , and X is again a diagonal matrices 
containing the eigenvalues of C. The good news is that 
these eigenvalues can be obtained through the application 
of the VEC[.] operator to the two-dimensional FFT of 
the autocovariance function c(ti,t 2 ) = E{J\f(j + T\,l + 
T 2 )7V(j, /)} of A/". Since also !F is a unitary matrix, it is 



x T C- l s » {x T T H ){TCT ti Y v {Ts 



H\-l, 



(B.5) 



= x H Y, 's = £"(DIAG[S "]0s), (B.6) 

Here, symbol "~" now indicates the two-dimensional FFT. 
Similar problems and solutions as in Sec. l2.ll h.old concern- 
ing the fact that, with the use of C, both S and X are 
implicitly assumed to be periodic functions along each di- 
mension with period N r and iV c , respectively. Once com- 
puted, filter u = C~ 1 s, that is a (N r N c ) x 1 array, can 
be converted in its original two-dimensional form simply 
columwise reordering its elements in a N r x iV c matrix U 
(i.e., the inverse operation of the VEC[.] operator). 



-H 



firstly by solving the system of linear equations = s, 
En ■ ■ ■ Sim 



E = 



(B.ll) 



'Nl ■ ■ ■ ±>MM 



that is highly structured and sparse (i.e., not computa- 
tional demanding), and then through x z. Alternatively, 
if x, s and n indicate the arrays (|B.7[) - (|B.9|) with the el- 
ements that are rearranged in row rollout order, and C_ 

is the covariance function of n, then T(x) = x S s, 
S » $C_$ H , can be efficiently computed through block- 
matrix operations exploiting the fact that S is a block 
diagonal matrix. 

Both approaches can be used to compute filters 
wsnr(^) and {7snr(^)- In particular, if the elements of 
x, s and n are arranged in row rollout order, a solution 
formally identical to Eq. (|A.12[) can be obtained if matrix 



(M) 



(B.12) 



is used in Eqs. (|A.8|) and (IA.11I). Again, this result co- 



inci des with that provided by iHerranz and Sanz 
and lHerranz et al. I (|2008bh . 



(2008a) 



B.2. Multiple-frequency observations 

In the case of multi-frequency observations, the situation 
becomes even worst since MF has to be applied to M 
signals at the same time. Again, through 

s = VEC [VEC[5i], VEC[5 2 ], . . . , VEC[5 M ]] ; (B.7) 
x = VEC [VEC[ATi], VEC[AT 2 ], . . . , VEC[Af A/ ]] ; (B.8) 
n = VEC [VEC[A/i], VEC[A/" 2 ], . . . ,VEC[Af M }] , (B.9) 

it is possible to obtain a problem that is formally identical 
to that given by Eq. ((2|). The only difference is that now 
C in Eq. fl27|) C is a (MN P ) x (MN p ) block matrix. In the 
case signals {Si} are two-dimensional iV r x N c maps, then 
each of the Cy blocks is constituted by a (N r N c ) x (N r N c ) 
BTTC matrix. In particular, Cu provides the autocovari- 
ance matrix of the ith image, whereas Cy, i ^ j, the 
cross-covariance matrix between the ith and the jth ones. 
If, again, each BTTB block Cy is approximated with a 
BCCB matrix Cy, then matrix 

5 = 7 (M) ®^, (B.10) 

can be used to diagonalize each of the blocks Cij in C. In 
this way, the statistics T(x) = x T C~ 1 s, can be computed 



